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ABSTRACT 

Standard formulations of smoothed particle hydrodynamics (SPH) are unable to re- 
solve mixing at fluid boundaries. We use an error and stability analysis of the gener- 
alised SPH equations of motion to prove that this is due to two distinct problems. The 
first is a leading order error in the momentum equation. This should decrease with 
increasing neighbour number, but does not because numerical instabilities cause the 
kernel to be irregularly sampled. We identify two important instabilities: the clumping 
instability and the banding instability, and we show that both are cured by a suitable 
choice of kernel. The second problem is the local mixing instability (LMI). This oc- 
curs as particles attempt to mix on the kernel scale, but are unable to due to entropy 
conservation. The result is a pressure discontinuity at boundaries that pushes fluids of 
different entropy apart. We cure the LMI by using a weighted density estimate that 
ensures that pressures are single valued throughout the flow. This also gives a better 
volume estimate for the particles, reducing errors in the continuity and momentum 
equations. We demonstrate mixing in our new Optimised Smoothed Particle Hydrody- 
namics (OSPH) scheme using a Kelvin Helmholtz instability (KHI) test with density 
contrast f :2, and the 'blob test' - a l:fO density ratio gas sphere in a wind tunnel - 
finding excellent agreement between OSPH and Eulerian codes. 

Key words: Multiphase Smoothed Particle Hydrodynamics, Numerical methods, 
Monte-Carlo methods 



1 INTRODUCTION 

Smoothed Particle Hydrodynamics (SPH) was first 
introduced as a tool for studying stellar structure 
l|Gingold fc MonaghanI 1 19771 : ILucvI Il977] 'l. but has since 
formd wide applic ation in all areas of theoretical astroph ysics 
l|Monaghanlll992ll. in engineering ([Libersky et al.ll 19931 ). and 
beyond (e.g. iHieber fc Koumoutsakosll200^ 

Although there are many varieties of SPH, the cen- 
tral idea is to repres ent a fluid by di s crete partic les that 
move with the flow (jMonaghanl Il992l : IPricd 120051 '). Typi- 
cally these particles represent the fluid exactly, though in 
some variants the fluid is advecte d on top of the particles 
l|Diltslll999l : lMaron fc H owes"2003''). The key advantages over 
Eulerian scheme£| are its Lagrangian nature that makes it 
Galilean invariant, and its particle nature that makes it easy 
to couple to the fast multipole method for gravity tha t scales 
as 0[N) (|Dehnenll2000l : iGreengard fc RokhlinHl987l ). How- 
ever, SPH has problems correctly integrating fluid instabil- 
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^ This does not apply to Lag ra ngian moving mesh schemes that 
are Galilean invariant llSpringe]|l2009l 'l . 



ities and mixing at boundaries (iMorrid 


1996bl: DiltsI 19991: 


iRitchie & Thomasll200ll: iMarri & Whitd 


2003': ' Agertz et all 



20011 ). Several different reasons have b een suggeste d for this 
in the literature so far. iMorrid (|l996bl ) and lPiltd (|lSj99l ) ar- 
gue that the problem owes to errors in the SPH gradients 
that do no t show good convergence for irregular particle dis- 
tributions. iPricd |2003) argue that the problem owes to the 
fact that entropies are discontinuous at boundaries, while 
the densities are smooth. This gives spurious pressure blips 
at boundaries that drive fluids of different entropy apart. 
They find that adding thermal conductivity at boundaries 
to smooth the entropies gives improved mixing in SPH. 
IWadsIev et al.l ( |2008 ) make a similar argument, phrasing the 
problem in terms of an inability for SPH particles to mix 
and generate entropy on the kernel scale. They find that 
adding a heat diffusion term to model subgrid turbulence 
gives improved mixing in SPH. Finally, Ritchie fc T homag 
(2001) suggest that the problem lies in the SPH density es- 
timate. They introduce a new temperature weighted density 
estimate that is designed to give smoother pressures at flow 
boundaries, thus combating the spurious boundary pressure 
blip. 

In this paper, we perform an error and stability analysis 
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of SPH in its most general form to understand wiiy mixing 
fails. In doing this, we show that all of the above authors cor- 
rectly identified one of two distinct problems with mixing in 
SPH. The firs t is an 0{h~^) e rror i n the momen tum equation 
identified by iMorrid (|l996bh and iDiltd l| 19991 ). The second 
relates to entropy conservation on t he kernel scale , as ad - 
dressed directly b vlPricd ||2007|) andlWadslev et all l|2008h . 
and indirectly by iRitchie fc Thornas ( 2001 ) . Having iden- 
tified the problem, we present a new method - Optimised 
Smoothed Particle Hydrodynamics (OSPH) - that, given 
sufficient resolution, correctly resolves multiphase fiuid fiow. 

This paper is organised as follows. In 321 a^nd 33 we 
briefly review standard SPH schemes and introduce our new 
OSPH scheme. We show that there are two distinct problems 
with mixing in SPH: the 'Eq error' in the momentum equa- 
tion, and the 'local mixing instability' (LMI), and we show 
how both can be cured. In ijH we present our impleme ntation 
of OSPH in the GASOLINE code (|Wadslev et al.ll2004l ). In fg] 
we use a Kelvin Helmholtz instability (KHI) test with den- 
sity contrast 1:2 and 1:8 to demonstrate mixing in OSPH. 
We show the effect of turning on each of the OSPH improve- 
ments one at a time, arriving at a solution that is in excel lent 
agreement with the Eulerian code RAMSES (|Tevssierll2002l ). In 
3ni we use the standard Sod shock tube test to demonstrate 
that OSPH can successfully model shocks. I n 33 w e revisit 
the 'blob test' introduced in lAgertz et al.l (|2007l ). finding 
excelle nt agreement betwe en OSPH and the Eulerian code 
FLASH (|Frvxell et all I2OO0I '). Finally, in 31] we present our 
conclusions. 



2 SMOOTHED PARTICLE HYDRODYNAMICS 

In SPH, the fiuid is represented by discrete particles that 
move with the fiow. The density of each particle is estimated 
by a weighted sum over its neighbours: 



JV 

j 



(1) 



where hi and rrij are the smoothing length and mass of parti- 
cle i and j, respectively; we define r^j = — Vj and similarly 
for other vectors; and 14^ is a symmetric kernel that obeys 
the normalisation condition: 



W{\r~r'\,h)d^r' 



and the property: 
lim W{\r-r'\,h) ^S{\r-r'\) 

h—^O 



In the limit N 



00, h 



(and using irij/pj 



(2) 



(3) 



d^r' 



equation [T] recovers the continuum flow density. 

The equations of motion for SPH are then derived by 
discretising the Euler equations - the continuity, momentum 
and energy equations: 

dt ^ 
dv _ VP 
dt p 

du P, 



dt 



-V-v 



(4) 
(5) 
(6) 



where p, v and u are the density, velocity and internal energy 
per unit mass of the flow, respectively. 

The Euler equations can be derived fro m the La- 
grangian for hydrodynamics (e.g. lBennettll2006l ): 



L = 



(7) 



and in many modern derivations of the equations of motion 
for SPH, equation[7]is discretised, rather than equations |4]|6l 
Replacing the volume ele ment dV wi th the volume per 
SPH particle m/p, we obtain (|Pricell2005l '): 



(8) 



and the standard SPH equations of motion then follow from 
the Euler-Lagrange equations: 
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(9) 



(10) 



(11) 



where Wij = W(|r,j|, /i,). 

Note that equation |9] is automatically satisfied by time 
derivative of the SPH density estimate (equation[T|. For this 
reason, equation [T] is often referred to as the integral form 
of the continuity equation. 

The above system of equations are closed by the equa- 
tion of state: 



(7 - 1) Pi 



(12) 



This standard approach to deriving the SPH equations 
of motion has the advantage that the resulting equations 
are coherent by construction - that is they are consistent 
with a Lagrangian. This gives very good conservation prop- 
erties for the fiow. It is also straightforward to calculate 
the necessary correction terms that arise if the smoothing 
lengths are a function of space arid tirn e h = h{r, t) (see e.g. 
iNelson fc Papaloizou|[l994l : |Pricell2005l ). We do not include 
these correction terms in this paper. 

However, this standard derivation leads to a scheme 
that cannot correctly model fluid mixing processes (see 31]), 
which motivates us to move to a more general derivation. 
Discretising each of the Euler equations separately leads to 
a free function for each equation: 77, (j> and C,, as well as a 
different smoothing kernel for each. This is the approach we 
take next in 33] In 333] and 333] we will then use an er- 
ror and stability analysis of these more general equations of 
motion to constrain the new functions t], (j) and C, and our 
new kernels. By choosing these new free functions and ker- 
nels such that they minimise the integration error, we will 
arrive at a new scheme that can, with sufficient resolution, 
correctly resolve multiphase fluid flow. 



Also called consistent llOeer et al.ll2007h . 
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3 OPTIMISED SMOOTHED PARTICLE 
HYDRODYNAMICS 

In the previous section, we presented a standard derivation 
of the SPH equations of motion. However, this standard 
derivation leads to a scheme that cannot correctly model 
fluid mixing processes (see iJT]). In this section, we move to a 
more general derivation of the SPH equations of motion. We 
show that, in general, we have a free function for each of the 
Euler equations: 77, <j) and as well as a different smoothing 
kernel for each. There is also a freedom in the energy equa- 
tion in the choice of integration variable (energy or entropy; 
tj3.2p . In tj3.3l and tj3.4l we will then use an error and sta- 
bility analysis of these more general equations of motion to 
constrain the new functions 77, <j) and C, and our new kernels. 
By choosing these new free functions and kernels such that 
they minimise the integration error, we will arrive at a new 
scheme that can, with sufficient resolution, correctly resolve 
multiphase fluid flow. 



3.1 A general derivation of SPH 

In general, we have some freedom in how we discretise the 
Euler equations (equations [4][6l) to obtain t h e equations 
of motion for SPH (see e.g. lMonagiianl 1 19921 : IPricd l2005l : 
iRosswodlioogl ). The gradients in the Euler equations can be 
expanded to include a new free function for each equation: 
T), (f) and (": 



dp 
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dt 
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(13) 
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In the continuum form, above, rj, (j) and C, cancel. But in the 
discrete SPH form, they remain giving a useful additional 
freedom dPricelboOsh: 
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(16) 
(17) 
(18) 



where Hij = [H{\rij\,hi) + H{\rij\, hj)] /2, Kij and Lij are 
symmetrised smoothing kernels - one for each Euler equa- 
tion. Standard SPH (SPH from here on) is a special case 



For this reason, we use instead a generalised integral form 
for the continuity equation: 



JV 

E 



which, taking the time derivative, gives: 
dpi 



N 

dt ^ 'rij 



where: 



(19) 



(20) 



(21) 



and7)=§. 

This reduces to the continuity equation (equation I16|) 
under the kernel constraint: HijVij = ViWij, and for e = 0. 
The latter can be satisfied by construction if rji — rjj (as is 
the case for SPH), or if 77 = 0. However, in the continuum 
limit (A'^ 00, /i — ^ 0), e — > and so e will vanish with in- 
creasing resolution. For this reason, eauation ll9l gives a valid 
approximation to the continuity equation for any choice of 
rj, with e simply contributing an additional error term. 



3.2 Energy versus entropy forms of SPH 

A final freedom in the equations motion for SPH comes 
from the energy equation. Equation [TS] is the standard 
energy form of SPH, but there is also an entropy form 
iGoodman fc Hernauistlll99ll : ISpringel fc Hernauistil20o3 ). 
Instead of the internal energy, u, we evolve a function A{s) 
- the entropy function - that is a monotonic function of the 
entropy s defined by the equation of state: 



P = Ms)p1 



(22) 



Away from shocks and in the absence of thermal sources 
or sinks, Ai is a constant of motion. Thus, taking the time 
derivative of equation [22] and substituting for equation 1121 
we recover: 
dui Pi dpi 



dt 



dt 



(23) 



by construction. Schemes that obey equation [23] are called 
thermodynamically consistent. 

In practice, we find - for the tests presented in this pa- 
per - that the energy and entropy forms of SPH give near- 
identical results, provided that equation [53] is satisfied (for 
adiabatic flow). We use the thermodynamically consistent 
energy form throughout this paper. This gives us the con- 



straints: C = ^ s-nd LijVi 



which we 



of the above with 77 



c 



apply from here on. We also use KijVij = ViWij, as in 
standard SPH. This is not a formal requirement, but en- 
sures that coherence is recovered in the limit of constant 
density. 



Equation [T5] casts the continuity equation in differential 
form. This is problematic since, in this case, the particles no 
longer represent the fluid exactly. Instead they represent a 
moving mesh on which the Euler equations are solved. This 
leads to the danger that high dens ity regions will contai n 
few particles leading to large errors aron fc Howes|[2003 ). 



3.3 Errors: choosing the free functions 

In this section, we perform an error analysis of the gener- 
alised equations for SPH (equations 1171 1181 and I19|l derived 
in ^3.11 We will then choose our free functions 77, (f) and so 
that these errors are minimised. 
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3.3.1 Error analysis 

We assume that the pressure and velocity of the flow are 
smooth. In this case, we can Taylor expand to give: 



(24) 



and 



V,- ~ Vi + /i(x,,- • V,)v, + 0(r) (25) 

where x^j = Vij/h, and we have assumed a constant smooth- 
ing length h. 

Substituting equations [24] and [25] into the continuity 
and momentum equations gives: 



^ ~ -p,; (R.VO ■ V, + e + 0{h) 
and 

5=.-AEo,.-(XiZi)^ + o(/.) 

at hpi pi 
where Eo,i is a dimensionless error vector given by: 



En 



JV 



(26) 
(27) 

(28) 



and Ri and are dimensionless error matrices given by: 

N N 

^ Pi 



V, = > S,, (29) 

^ Pi 



with: 



^ _ 1 dW,, 

^ii — 



X dx 



2 

yijXij yijZij 



(30) 



Pj Vi 



where Vf = ftVi; Xij = {xij,y,j, z,j)- x = |xij|; /ij = 
and (7i, — 

The accuracy of the continuity equation (|26|) is given by 
the extent to which e = (see equation I2ip and Ri = I, the 
identity matrix. The accuracy of the momentum equation 
l|27p is given by the extent to which Eo,i = and Vi = I. 
(The energy equation behaves similarly to the continuity 
equation with e = 0.) 



/ = 



P(r') 
P(r) 



1 + /i 



(x- 



Vp + 0{h^ 



(33) 



and we see that, by symmetry of W, R = I to 0{h?). In 
fact, Taylor expanding to an order higher than above, it is 
straightforward to show that the whole continuity eq uation 
is ac curate to 0{h?) in the limit N oo (see e.g. jPricel 
I2OO5I ). A similar argument applies to the other SPH equa- 
tions of motion and leads to the conclusion that SPH is 
accurate to 0{h?). However - and this is a key point - this 
formal calculation is only valid for smoothly distributed par- 
ticles tn the limit N 00. \n practical situations, where 
we have a finite number of particles within the kernel and 
these are not perfectly smoothly distributed, the leading or- 
der errors in the continuity equation appear at 0(0) and are 
contained within the matrix R. We will quote orders of error 
from here on in this finite particle limit. 

We can think of each term of R as a finite sum approx- 
imation to a dimensionless integral that should be either 
(for the off diagonal terms), or 1 (for the diagonal terms). 
For smooth particle distributions, this approximation is a 
good one since fij ~ 1, while mjj pj gives a good estimate 
of the volume of each particle within the kernel. However, if 
the particles are distributed irregularly on the kernel scale - 
for example at a sharp density step - then /ij can grow arbi- 
trarily large, while ruj/ pj becomes a poor volume estimate. 
We will demonstrate this in 33 

We can improve matters by choosing rj = p, which fixes 
f — 1 always. However, the integral form of the continuity 
equation then becomes: 



JV 

pi = —mjWij 



(34) 



which must be solved iteratively and is not guaranteed to 
converge. Worse still, e is now no longer zero and contributes 

an additional error. 

[Ritchie fc Thornai (|200l[ ) present an interesting solu- 
tion to this dilemma. If the pressures are approximately con- 
stant across the kernel (Pi ~ Pj) then, for the energy form 
of SPH (see equation [T^ and >i3.2p . ^ ~ and equation 
1341 is well approximated by the integral continuity equation: 



3.3.2 Minimising errors: the continuity equation 

Let us consider how accurately equation [26] approximates 
its Euler equation equivalent (equation [3]| . First, consider 
standard SPH where r] = 1 and e = by construction. 
Typically in the literature, the error is calc ulated only in 
the continuum limit {N — >■ oo\ ft — > 0; see e.g. |Pric3 l|2005h ). 
In this case, the sums become integrals, and (using ruj/ pj — >■ 
d?x') we obtain terms like: 



d^^7(x,x')ii^f: 
X — x' ox 



(31) 



lim R33(x) : 
and 

hm Ri2(x) = / da;/(x,x) ■ — (32) 

JV^oo Jv |X — X'l OX 

where the notation 33 refers to element [3, 3] in the matrix 
R. 

If we assume smooth densities, then we can Taylor ex- 
pand / also to obtain: 



JV 

Pi = > —rUjWi 
^ Ui ' 

j 



(35) 



This can be solved without the need for iteration. 

The above suggests that we use r) — 1/u. Thermody- 
namic consistency then requires that we set — rj — 1/u 
(see i32}. 

There may be some advantage, however, to using the 
entropy form of SPH. In this case, the equation of state is 
given by equation 1221 For approximately constant pressure 



across the kernel, we now have that — 

' Pj 

integral continuity equation becomes: 

'A, 



Pi 



JV 

E 



nijWij 



and the 



(36) 



This has the advantage that, in the absence of shocks or 
thermal sources/sinks, ^ = and so the error term e = by 
construction (see equationl2ip. In practice, however, we find 
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no appreciable difference between the energy and entropy 
forms of SPH for the tests presented in this paper. This 
suggests that e is not a significant source of error. 

Equations 1351 and 1361 retain the desirable integral form 
for the density, while giving significantly improved error 
properties. They also have a second important advantage 
that we discuss in i|3.5l We refer to equation [35] as the 'RT' 
density estimator for the energy form of SPH; and equation 
[36l as the RT density estimator in the entropy form. 

3.3.3 Minimising errors: the momentum equation 

The momentum equation (eauation l27|) is more problematic 
than the continuity equation. Its accuracy is governed not 
only by the extent to which Vi — I, but primarily by the 
leading Eo,i term that should vanish. 

First, consider the situation in standard SPH where 
(f) — 1 and gij = Pj/pi. As for the continuity equation, in 
the continuum limit (TV — > oo;h — > 0), Eo = + O{h^) since 
\/fWij is antisymmetric. However, this analysis is only rel- 
evant if the particles are smoothly distributed on the kernel 
scale. For irregularly distributed particles, gij can grow ar- 
bitrarily large, while rrij/pj is not guaranteed to be a good 
volume estimate. In such situations, Eo contributes a signif- 
icant error. Worse still, moving to higher resolution is not 
guaranteed to help. In order for the SPH integration to con- 
verge as ft 0, we require that Eo,i shrinks faster than 
. This requires some care in making sure that h does not 
shrink too fast as the number of particles is increased. 

A density step is an extreme example of an irregular 
particle distribution, and this suggests that the Eq error is 
at least in part responsible for SPH's failure to correctly 
model mixing processes between different fluid phases. We 
demonstrate this in !j5] 

There are three key problems with ensuring that 
j0-Eo,i will shrink with increasing resolution. The first is 
the function gij. In SPH, this is the ratio pj/pi which is 
large when there are large density gradients. We can signif- 
icantly improve on this we choose our free function (j) = p. 
In this case, we have gij — g~.^ = 1 by construction, and gij 
no longer contributes to the Eo error even for large density 
gradients across the kernelf]. The second problem relates to 
kernel scale smoothness. If particles clump or band on the 
kernel scale, then we will have poor kernel sampling and Eo 
will not approach its integral limit even at very high reso- 
lution. Ensuring that this does not happen means ensuring 
that our OSPH scheme is stable to perturbations. We dis- 
cuss this next in tj3.4l The third and final problem is the 
volume estimate of each particle mj/ pj. This will be poor if 
the particles are irregularly distributed on the kernel scale 
(for example at a density step) leading to a large Eo error. 
We discuss this further in H3.5I 

^ It is interesting to note that other work in the literature has 
also foun d that tj!) = p is the preferred choice if density gradients 
are l arge jOger et al.ll2007l : |Pricj|2005l : iRitchie fc ThomajlioOll : 
iDiitsHigggT But we could not find a detailed proof similar to that 
presented here. Interestingly, iMarri & White! ||2003| ) find empir- 
ically that = p'^/^ gives the best performance for multiphase 
flow in their tests. Our analytic results here suggest that this is 
not the optimal choice, though perhaps the inclusion of cooling 
and/or other physics makes a difference. 



The choices — rj = 1/u; (p — p and the kernel con- 
straints HijYij = LijYij — KijVij — ViWij are the first 
important ingredients in our OSPH scheme. These choices 
mean that we are no longer coheren10, but this only in- 
troduces tolerable O(h^) errors in the energy conservation 
(|Hernauist fc KatjflQSgT ). 

3.4 Stability: the choice of kernel function 

In ^3.31 we used an error analysis of the generalised SPH 
equations of motion to show that the dominant source of 
error in SPH is in the momentum equation - the Eo error. 
We showed that choosing the free functions C = ^ = 
<j} ~ p and the kernel constraints HijVij — LijVij = KijVij = 
ViWij should minimise both this error and errors in the 
continuity equation, and we called these choices OSPH. 

In OSPH, provided the particles are regularly dis- 
tributed on the kernel scale, we can make Eo arbitrarily 
small simply by increasing the neighbour number. However, 
if the particles are irregularly distributed, Eo can shrink 
very slowly with increasing resolution. In this section, we 
show that for large neighbour number the cubic spline ker- 
nel typically used in SPH calculations is unstable to both 
particle clumping ( H3.4.ip and particle banding (|3.4.2|l . and 
we derive a new class of kernels that are stable to both even 
for large neighbour number. In !j5] we show that these new 
kernels give significantly improved performance at no addi- 
tional computational cost. 

3.4.1 The clumping instability 

The clumping instabilitjjf] can be derived from a linear 3D 
sta bility analysis of t he OSPH equa tions of motion. Follow- 
ing [Morri^ [19963 and lMorrislll996bl . we imagine a lattice of 
equal masses m of equal separation (Aa;o, Aj/o, Azo) with 
initial density po and pressure Po- We perturb these with a 
linear wave of the form: 

Xi = xo,i -I- aexp[i(k ■ xo,i — (37) 

pi = P0 + D exp[i(k • xo,i - a;t)] (38) 

Pr = Po + cId exp[i(k ■ xo,» - ojt)] (39) 

and similar for particle j, where = ^ = 'yPo/po is the 
sound speed assuming an adiabatic equation of state (7 = 1 
gives an isothermal equation of state), a = {X,Y,Z) is the 
amplitude of the perturbation and k = {kx,ky,kz) is the 
wave vector. 

To simplify the analysis, we assume that we have a 
lattice symmetry such that for every displacement vec- 
tor xo,ij — xo.i — xo,j to a neighbour, there is also one 

* Note that it is possible to construct pseudo-coherent versions 
of OSPH using ( = ri = (p = 1/u for the energy form, or 

f = »7 = = for the entropy form . Introducing 'grad 

h' terms as in iNelson fc Papaloizoul lll994l '). such schemes can 
then be made to conserve energy exactly in the limit of constant 
timesteps. However, they are only truly coherent up to the ap- 
proximation that the e error in the continuity equation is small 
(see equation [2TJ . Nonetheless, it would be interesting to explore 
such schemes in future work. 
^ Also called the tensile instability. 
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Cubic Spline (CS) kernel: 
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Figure 1. Stability plots for the Cubic Spline (CS) kernel (equation 144 l l in OSPH. The plots show contours of the frequency a;^/fc^/c^ 
of plane waves impacting a regular lattice of particles, as a function of the wavenumber k and the smoothing length /i, in units of the 
inter-particle spacing dx = 1. From left to right the plots show {kx,ky,k^) = fc(l, 0, 0), fc(l, 1, 0) and fc(l,l,l). The three rows show 
the longitudinal wave and the two transverse waves for each of these orientations. Also marked by the blue dashed lines are the h that 
corresponds to 32 neighbours (bottom line) and 128 neighbours (top line). OSPH is unstable if < (grey regions). 



at — xo.ij. Then, plugging equations 1371 1381 and [39] in to 
1171 discarding terms higher than first order, and connect- 
ing D to X, Y, Z through the continuity equatioijf] [D = 



dx^ dy 

the 3D OSPH dispersion relatiorQ: 



), we obtain 



° We use here the full OSPH continuity equation in the entropy 
form (equation 1161 with rj = 1/A^i'''). However, for plane waves 
on a constant density lattice Aj/Ai = 1 and so this is identical 
to the SPH continuity equation with r] = 1. 



^ This is actually identical to the S PH dispersion r elation derived 
under the same assumptions in 3D llMorrislll996lJ) . 
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Core Triangle (CT) kernel: 
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Figure 2. As Figure [T] but for the Core Triangle (CT) kernel. 



2mPo 



^H(Wo,»j)(l -cosk-xo,ii)+ 



(7-2)^(q,AqO 
Po 



where A qi is the outer product of q^, H(W^) is the Hes- 
siailfl of W: 



Haa 
Hab 



d^W _ d'W(r) ^ dWjr) 1 ^ ^ Xa 



dXa^ 

dW 



dr r 



d^Wir) Xaxt dW{r) XaXl) 



dxadxb dr^ 
and qi is given by: 



dr 



(41) 



(42) 



(40) 



° Recall that the outer product of two vectors is a matrix, while 
the Hessian is a square matrix of second order partial derivatives. 
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High Order Core Triangle HOCT4 kernel: 




TS — ^ — I — 5 — ! — 
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-0.01 
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Figure 3. As Figure[T] but for the High Order Core Triangle (HOCT) kernel with n^; = 4 



T—\ — 5 — 5 — I — s~~r 



Qi = ^sink • xo,ijViVKo,ii (43) 
j 

Our scheme is stable if oj'^ ^ 0. It is also desirable for 
the numerical sound speed to equal the true sound speed: 

2 II 2 2 

In SPH it is typical to use the cubic spline (CS) kernel 
given by: 

( l-&x^ +&x^ < a; i 
W = 2(l-xf i < a; < 1 (44) 

[ otherwise 



where x = r/h is the distance from the centre of the kernel 
in units of the smoothing length. 

In Figure[T] we show contours of cj^/fc^/c^ as a function 
of wavenumber k and the smoothing length h in units of the 
inter-particle spacing dx — 1, for the CS kernel. We assume 
an adiabatic equation of state with 7 — 5/3. From left to 
right the plots show {kx,ky,kz) = fc(l, 0, 0), fc(l, 1, 0) and 
fc(l, 1, 1). The three rows show the longitudinal wave and 
the two transverse waves for these orientations. 

From Figure [T] it is clear that the CS kernel in 3D is 
unstable to longitudinal waves for h ^ 2, and very unsta- 
ble to transverse waves. The unstable longitudinal waves 
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drive the clumping, or tensile i nstability that causes parti- 
cles to clump on the ke r nel sc a le jSchuessler fc Sclimittll981 
Thomas fc CouchmanI Il992l : iHerantI 1 19941 : I Morris! Il996a 
Monaghan.2000. ). 



The clumping instability is a problem because it means 
that increasing the neighbour number will not give improved 
sampling of the kernel, and the Eq error will remain large. 
However, the situation is dramatically improved if we add a 
constant central core to the kernel gradient This gives 
a constant force term at the centre of the Kernel that phys- 
ically prevents clumping. We choose a kernel that is maxi- 
mally similar to the CS kernel, while obeying — const. 
Vr < a, where a is the core size. This leads us to the Core 
Triangle (CT) kernel: 



W = 



iV 

ft3 



(-12a + 18a'^) X + P < 2: < a 

1 - 6a;^ + 6a;^ a < a; | 

2(1 -a;)^ i<a;<l 

otherwise 



(45) 



where /3 = 1 + 6q^ - 12Q^ iV = 8/[n (6.4a '"^ - 16a® + l)], 
and the core size is fixed at a = 1/3 by the requirement that 
be continuous. 

Figure [2] shows stability plots for the CT kernel. The 
CT kernel has greatly improved stability for the longitudi- 
nal waves (top row) compared to the CS kernel and should 
give significantly improved performance for large neighbour 
number. We demonstrate this in iJS] 

Note that for all of the kernels we use in this paper, we 
consistently apply the kernel for the density estimate and its 
gradient for the energy and momentum equations. For small 
neighbour number, the central triangle in the CT kernel will 
degrade the quality of the density estimate. However, in this 
paper we typically use large neighbour numbers (> 100). In 
this case, very few particles sample the inner regions of the 
kernel and the bias introduced in the density is negligible. 
(The quality of the density estimate in OSPH can be seen 
in the Sod shock tube test in ^) We found in tests that 
retaining the CS kernel just for the density estimate gives 
near-identical results. 



3.4-2 The banding instability 

The clumping instability is a result of unstable longitudinal 
waves. A related instability - the banding instability - is a 
result of unstable transverse waves. For both the CT and CS 
kernels, there are broad bands of instability to transverse 
waves (see Figures [1] and [2]). If the neighbour number is 
carefully chosen to lie in a stable region, banding will not 
occur. However, banding can still be excited at boundaries 
if h changes there, moving into an unstable region. 

For both the CS and CT kernels, it is difficult to find a 
suitable neighbour number for which the kernel is stable to 
all transverse and longitudinal modes. This suggests hunting 



for an even more stable kernel. A full search is beyond the 
scope of this paper. Here, we present a simple class of kernels 
that improve sta bility by moving to higher order (|Morrisl 
[l996b. V Following IPricel C2005), we generalise our CT kernel 
to order Uk to obtain the following class of kernels that we 
call the High Order Core- Triangle (HOCT) kernels: 



Px + Q 
(1 - a;)"*- 
B{P - x) 
{l-xY" 

{i~xY<^ 





+ ^(7- 

+ A(7- 



< X ^ a 
a < X ^ j3 

/? < X < 7 
7 < a; ^ 1 
otherwise 



(46) 



where: 



A = 



1 



B 



7"fc-3('y2 _ ^2) 

1-f- A7"'=-i 



(47) 



(48) 



P = -nfe(l-a)"''-'-nfcA(7-a)"'=-'-nfcB(/3-a)"'=-'(49) 

Q={1- aT" + ^(7 - a)"" + B(/? - a)"" - Pa (50) 

and a and TV are calculated numerically for a given choice 
of nfe. Continuity requires that a solves the equation: 



= (1 - a) 



-A(7- 



+ B(/3- 



(51) 



where /? and 7 are free parameters. In this paper we choose 
P = 0.5, 7 = 0.75. Other choices, and indeed other high- 
order kernels, may give better results than those presented 
here. We tabulate values for A, B, P, Q, a and A'^ as a func- 
tion rifc in Table[T] Notice that the core size a decreases with 
rik- 

Stability plots for the HOCT4 kernel (with Uk = 4) 
are given in Figure [S] Notice the improvement over the CT 
kernel, particularly for the transverse waves. There are two 
bands where the kernel is fully stable to both longitudinal 
and transverse waves on a lattice: 96 neighbours and 442 
neighbours, corresponding to h — 2.86 and h = 4.75, re- 
spectively. We use the latter choice since this also gives very 
low Eq. The CT kernel also has a stability band for h ~ 4.75, 
but this is narrower than for the HOCT4 kernel, while the 
H0CT4 kernel with this many neighbours gives better spa- 
tial resolution. (It is important to realise that the smoothing 
length for different kernels takes on a different meaning in 
terms of spatial resolution. We suggest a resolution criteria 
based on the numerical sound speed versus the true sound 
speed for longitudinal waves: ai^/fe^/c^. Spatial scales are 
well resolved if uj^ /k'^ /cg ~ 1. By this definition, our choice 
of 442 neighbours [h = 4.75) for the H0CT4 kernel gives a 
very similar spatial resolution to 128 neighbours {h = 3.14) 
for the CT kernel.) 

Note that our stability analysis only applies for particles 
arranged on a lattice. Hexagonal close-packed particles, ran- 
domly arranged particles, and indeed boundaries may have 
different preferred stability regions. A full analysis is beyond 
the scope of this present work. 

The banding instability is not as problematic as the 
clumping instability for the tests we present in ^ and ^ 
Unlike the clumping instability, it does not seem to (directly) 
play a major role in preventing mixing from occurring in 
SPH (see !j522}. 
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SPH-CS-128 




Figure 4. A Kelvin-Helmholtz instability (density ratio Rp = 2) at trh = 1 modelled with SPH, TSPH and OSPH using CS, CT and 
HOCT4 kernels (see equations 1441 l45l and l46ll . From left to right the plots show, in a slice of width dx = 1 about the z-axis: density 
contours; a zoom-in on the particle distribution around one of the rolls; the magnitude of the |Eo| error (see equation I28|l as a function 
of y \ and the pressure in a slice of width dx = 1 about the x-axis, as a function of y. The circles on the density contour plots mark the 
size of the smoothing kernel, h. 
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Figure 5. A zoom in on the boundary for the KH test in (from left to right): TSPH-CT-128, TSPH-HOCT4-442 and OSPH-HOCT4-442 
S't tkh = 1- The plots are vertical projections of all the points within a 64 X 64 X 2 cuboid. Each point is plotted as a solid black circle 
with some transparency in order for the gaps to stand out. 



Table 1. Parameters for a selection of High Order Core-Triangle 
(HOCT) kernels. See equation 1461 for details and definitions. 





A 


B 


P 


Q 


a 


N 


3 


2.4 


-9.4 


-1.81 


1.028 


0.317 


3.71 


4 


3.2 


-18.8 


-2.15 


0.98 


0.214 


6.52 


5 


4.27 


-37.6 


-2.56 


0.962 


0.161 


10.4 


8 


10.1 


-300.8 


-3.86 


0.942 


0.0927 


30.75 



3.5 The local mixing instability & RT densities 



Our error analysis in ii3.3l missed one very important error 
term. This is because the Taylor expansion assumed that 
both the pressures and velocities in the flow are smooth. 
Unfortunately, in SPH at sharp boundary this is not the 
case. The reason for this is easiest to understand using the 
entropy form of SPH, as follows (similar arguments apply 
also for the energy form). 

Imagine a density step of ratio Rp = pxj pi initially in 
pressure equilibrium, such that the entropy function (equa- 
tion [221) is given by AxjAi = l/R],. Now imagine that we 
perturb the boundary very slightly by pushing a low density 
particle towards it. The particle's entropy is conserved, but 
its density increases very rapidly proportional to Rp. This 
leads to an increase in pressure: Pi = P -I- k,\R, where ki is 
some constant that depends on the perturbation size and the 
kernel. On the other side of the boundary, if we push a high 
density particle towards the low density region, however, its 
density will rapidly decrease giving a decrease in pressure: 
P2 = P — K2Rp. This drives us towards a pressure discon- 
tinuity at the boundary which drives an associated error in 
the momentum equation. It can be thought of as a funda- 
mental result of particles trying to mix on the kernel scale, 
but being unable to as a result of entropy conservation. We 
call this the local mixing instability (LMI). 

Although not phrased in terms of the LMI, the LMI 
is a recognised problem in the literature and there are es- 
sentially two classes of solution. We can generate entropy at 
the boundary to give smooth entropies and therefore smo oth 
pressures, as in jWadslev et al.ll2008l ') and iPricd (|2007l ): or 
we can try to obtain sharper densities that are consistent 
with the discrete entropie s. This is the approach adopted by 
I Ritchie fc Thomad (|200ll '). and the approach we take in this 



paper. The key advantage of sharpening the densities is that 
we do not need to specify a subgrid mixing model. 

The sharper densities we require are exactly what we get 
from the density estimate given in equatio n 1361 and origi- 
nally proposed bv lRitchie fc Thomad |200ll ) - the 'RT' den- 
sity estimate. Consider the perturbation discussed above, 
but now using the RT density estimatE0. A low density par- 
ticle which has half of its kernel in the high density phase 
(an extreme example), will have a density: 



3 



my 



(52) 



where A'^; is the number of particles in the low density region, 
Nr is the number in the high density region, and we have 
used the fact that the ratio (Aj/Ai) 1 = 1 for the low density 
region. 

If the simulation is adiabatic and started in pressure 
equilibrium, then for the high density region [Aj/Ai)~ — 
1/Rp, and since the high density particles sample the ker- 
nel Rp times more often than the low density particles, we 
recover: 



Plo 



(53) 



which is identical to a particle in the low density region. 
A similar derivation applies to a high density particle at 
the boundary. Thus, the RT density estimate ensures that 
densities remain sharp. It is straightforward to show that it 
also ensures the pressures are single valued throughout the 
flow. Substituting the RT density estimator (equation I36|) 
into equation 1221 we obtain: 



■ JV 
. i 



(54) 



Notice that the entropy function Aj now appears inside 
the sum, whereas in standard SPH it would appear as Ai 
outside of the sum. This difference ensures that the Pi will 
be everywhere single valued throughout the flow - even at 

^ We use here the entropy form given in equation I36l since we use 
the entropy form of SPH in this analysis. If instead, we use the 
energy form of SPH then we should use instead equation 1351 
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TSPH-CT-128 




OSPH-HOCT4-442; 1:8 density ratio 




Figure 6. Long term evolution of the KH instability in TSPH and OSPH versus the Eulerian code RAMSES. Prom left to right, the panels 
show density contours in a slice of width dx = 1 about the z-axis at times tkh = li 2 and 3. 
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Figure 7. A Sod shock tube test in SPH (top) and OSPH (bottom). From left to right the panels show the variation in density, velocity, 
pressure and temperature across the shock, respectively. The blue lines give the analytic solution. This test was performed in 3D. 



boundaries. A similar derivation can be made for the en- 
ergy form of SPH, in which case we should use the density 
estimate given in equation 1351 

The RT density estimate is robust to particle mixing on 
the kernel scale and should lead to a dramatically reduced 
LMI. We demonstrate this in ij5l Furthermore, the RT den- 
sity estimate ensures that our error analysis in i]3.3l is valid 
by construction since it ensures smooth pressures (recall that 
we assumed that both the pressures and the velocities were 
smooth, but not the densities). And, since the RT density 
estimate leads to sharper densities, it gives improved vol- 
ume estimates for the particles. This suggests that we can 
expect the RT density estimate to reduce Eq at boundaries. 
We demonstrate this also in 33 

Note that the RT density estimate is chosen to ensure 
single valued pressures throughout the flow. However, when 
extracting results from a simulation, it is the positions of 
the particles themselves that describe the state of the fluid. 
This suggests using the density estimate in equation [1] for 
calculating the observable flow density, rather than the RT 
density estimate. This is the approach we adopt in this pa- 
per, though the difference is negligible. 



4 IMPLEMENTATION 

We implemen ted ^OSPH in the GASOLINE code 

l|Wadslev et al.l |2004 ). a parallel implementation of 
TreeSPH that uses a fixed number A'^ of smoothing neigh- 
bourf^. and a standard prescriptio n for the artificial 
viscosity as in iGingold fc Monaghanl (|l983D with Q = 1, 

A llowing for varying neighbour number is needles sly dissipa- 
tive llNelson fc PapaIoi2oulll994 lAttwood et al.ll2007t) . 



(3 = 2, controlled with a Balsara switch l|Balsaral Il989l ). 
We used variable timesteps controlled by the Courant time 
with a Courant factor of 0.4. 

The improved stability and error properties of OSPH 
motivate a full re-examination of the standard SPH artificial 
viscosity. This is beyond the scope of this present work. How- 
ever, we note that the improved stability in OSPH means 
that particles better follow characteristics of the flow, while 
the gradients in the Balsara switch will be less noisy. Both of 
these effects should act to decrease the viscosity in regions of 
steady flow. (Note that all numerical schemes carry numer- 
ical viscosity, whether it is manifested through limited res- 
olution or artiflcial shock-capturing viscosity. Indeed, these 
viscous terms are vital for successfully modelling shocks.) In 
35] ^and ^ we show that our OSPH results agree very well 
with analytic expectations, and with the results from Eule- 
rian codes. This suggests that the viscosity prescription in 
OSPH is not a signiflcant source of error. Certainly it is not 
responsible for SPH's inability to model mixing processes. 



5 THE KELVIN-HELMHOLTZ INSTABILITY 

In this section, we use a 1:2 and 1:8 density ratio shearing 
fluid simulation to test mixing in OSPH. We use the naming 
convention XSPH-K-N, where X denotes the variety of SPH, 
K the choice of kernel, and N the neighbour number (see 
Table [2]). 



5.1 Numerical set-up 

A Kelvin-Helmholtz instability (KHI) occurs when two 
shearing fluids are subjected to an infinitesimal perturba- 
tion at the boundary layer. The result of the perturbation 
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Figure 8. The blob test in SPH, TSPH, OSPH and the Eulerian code FLASH. From left to right the plots show density contours at times 
TKH = 0, 1, 2 and 3. The contour bar gives logarithmic density in cgi units. 



is a linearly growing phase in which the layers start to in- 
terpenetrate each other, progressively developing into a vor- 
tex in the non-linear phase that mixes the two fluid layers. 
The growth-rate of the instability is in general a complicated 
function of the shear velocity, fluid densities, compressibility, 
interface thickness, gravity, viscosity, surface tension, mag- 
netic fleld strength etc. In this test, we are only interested 
in the behaviour of inviscid, incompressible (i.e. with bulk 
motions very much less than the sound speed) perfect fluids 



neglecti ng gravity. In this cas e, the linear growth rate of the 
KHI is l|Chandrasekhaj|l96ll ): 



w = k 



(piP2 



,1/2, 



(Pi + P2) 



(55) 



where fc = 27r/A is the wavenumber of the instability, pi and 
P2 are the densities of the respective layers and v = vi ~ V2 
is the relative shear velocity. The characteristic growth time 
for the KHI is then: 
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Figure 9. The magnitude of the |Eo | error in a slice of width dx = 1 about the z-axis, as a function of y (top); and the pressure in a sHce 
of width dx about the x-axis and dx about the z-axis, as a function of y (bottom) for the blob test in SPH-32 (loft), TSPH-HOCT4-442 
(middle) and OSPH-HOCT4-442 (right) at time tkh = 1- 



Flavour of SPH 



c 



Kernel 



SPH 111 CS 

TSPH 1 p 1 CS, CT, HOCT4 

OSPH 1/u p 1/u HOCT4 

Table 2. The different flavours of SPH we explore in this work. 
The free functions rf, <j> and C, are defined in equations [161 [TT] and 
1181 The Cubic Spline (CS) kernel is given by equation 1441 the 
Core- Triangle (CT) kernel is given by eauation l45l and the fourth 
order High-Oder Core Triangle (HOCT4) kernel is given by equa- 
tion [46] with rifc = 4. 



T"KH 



Stt _ (pi + P2)A 



(56) 



This is a particularly challenging test for SPH/OSPH be- 
cause the velocity due to particle noise can approach the 
sound speed which can wash out the physical velocity per- 
turbation relevant for this test. 

For the simulations, we set up the problem in three 
dimensions using a periodic thin slab defined hy x £ 
{-0.5,0.5}, y e {-0.5,0.5} and z e {-1/64,1/64}. The 
domain satisfied: 



pi,Ti,vi 

P2,T2, V2 



\y\ < 0.25 
\y\ > 0.25 



(57) 



The density and temperature ratio were Rp — pi/p2 = 
T2/T1 — C2/C1, ensuring that the whole system was pressure 
equilibrium. The two layers were given constant and oppos- 
ing shearing velocities, with the low density layer moving at 
a Mach number Ai2 ~ —V2/C2 ~ 0.11 and the dense layer 
moving at AI1 = Rp. The density ratios considered in 

this work are small which assures a subsonic regime where 
the growth of inst abilities can be treated using equation 1561 
(jVietri et aLlllQQTh . 

To trigger instabilities, velocity perturbations were im- 
posed on the two boundaries of the form: 



5uy[sin(27r(j: -f A/2)/A) exp(-(10(y - 
- sin(27r2:/A) exp(-(10(y 4- 0.25))^)] 



0.25))^ 



(58) 



where the perturbation velocity 5vy /v — 1/8 and A = 0.5 is 
the wavelength of the mode. 

Equal mass particles were placed in lattice configura- 
tions to satisfy the setup described above. To satisfy pressure 
equilibrium everywhere, in TSPH the temperatures were ad- 
justed at boundaries to be coherent with the smoothed den- 
sity step measured by equation [T] This was not done for the 
OSPH simulations since these sharpen the densities using 
the discrete initial temperatures. 

The low density region p2 was set up using 256 particles 
in the a;-direction and the appropriate number of particles 
in the other dimensions to satisfy a fixed inter-particle dis- 
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tance. The high density region pi was created in the same 
way with 320 particles in the x-direction. We adopted a pe- 
riodic simulation domain. 

The RAMSES simulation used the same numerical set-up 
as described above, but in 2D rather than in a thin slab. We 
perfor med the Rp — 2 simulation using the LLF Riemann 
solver (|Tordl 19991 ) on a 256 x 256 fixed Cartesian grid. The 
LLF solver is rather diffusive and is used in order to suppress 
the growth of undesirable small scale KHIs arising from grid 
irregularities. 

We note that all numerical schemes carry numerical vis- 
cosity, whether it is manifested through limited resolution or 
artificial shock-capturing viscosity. A detailed study of this 
effect on the KHl and the relation to physical viscosity is 
beyond the scope of this paper. 

5.2 Results 

Figure [4] shows our results for the KHI test (density ratio 
Rp = 2) at TKH = 1 modelled with SPH, TSPH and OSPH, 
using three different kernels: CS, CT and HOCT4, and dif- 
ferent neighbour numbers as marked on each plot (see also 
Table [2]). From left to right, the panels show, in a slice of 
width dx — 1 about the z-axis: density contours of the simu- 
lation box, a zoom in on the particle distribution around on 
of the rolls, the magnitude of the |Eo| error (equation [28]) 
as a function of y, and the pressure as a function of j/ in a 
slice of width dx = 1 about the x-axis. 



5.2.1 The clumping instability 

Using the standard CS kernel, SPIf-CS-128 (top row. Fig- 
ure [4| and TSPH-CS-128 (second row) give poor results 
that improve very slowly with increasing neighbour num- 
ber. This can be seen both in the lack of strong evolution 
on the boundary, and in the large |Eo| error, even for 128 
neighbours. TSPH-CS-128 gives slightly better results than 
SPH-CS-128, showing the first beginnings of a KHI roll, but 
both are in poor agreement with the RAMSES results (bottom 
row). 

The reason for the poor performance in both SPH-CS- 
128 and TSPH-CS-128 is the clumping instability ( ^3XT|) . 
Particles gather together on the kernel scale, giving poor 
kernel sampling, and poor associated error. This can be seen 
in the particle distribution for SPH-CS-128 and TSPH-CS- 
128 (second row, Figure[4| which show visible holes and over- 
densities in the particle distribution. Using instead the CT 
kernel introduced in i]3.4.1l the results improve dramatically 
(third row. Figure |4]). Now the errors reduce for increasing 
neighbour number (see Appendix |A]l- With 128 neighbours, 
we successfully resolve a KH roll up to tkh = 1 with the 
correct growth time. 

It has been noted previously in the literature that 
putting a small core inside a cu bic spline kernel sup- 
press e s the clump ing instability (|Thomas &: CouchmanI 
I1992I : iHerand Il994l ). though its importance for modelling 
multiphase flow was not realised. Alternative fixes include 
adding an negative pressure term iMonagha nl l200d ). which 
in tests we find works also. However, we prefer changing 
the kernel to introducing new forces since we may then still 
estimate our errors through lEo|. 



5.2.2 The banding instability 

In addition to the clumping instability, there is also an insta- 
bility to transverse waves - the banding instability ( i]3.4.2p . 
For the KHI tests we present here, the banding instability 
occurs only on the boundary and appears to be relatively 
benign. This is shown in Figure [5] that shows a zoom in on 
the boundary at tkh = 1 for TSPH-CT-128, TSPH-H0CT4- 
442 and OSPH-HOCT4-442. The TSPH-CT-128 simulation 
has a kernel and neighbour number combination that are 
unstable to transverse waves (see Figure [2]), and banding is 
clearly visible on the boundary. However, TSPH-H0CT4- 
442 should be stable to transverse waves, yet the banding 
persists. Only in our full scheme, OSPH-HOCT4-442, is the 
banding is gone. 

To understand the above results, we ran an additional 
test that we omit for brevity - TSPH-HOCT4-96. This sim- 
ulation showed little boundary evolution because the low 
neighbour number and associated large |Eo| significantly 
damped the KHI. However, interestingly, there was no band- 
ing observed on the boundary (recall that 96 neighbours for 
the HOCT4 kernel should be stable to both transverse and 
longitudinal wave perturbations). 

Taken together, our results suggest that the observed 
banding at the boundary is a result of a transverse wave in- 
stability driven by the local mixing instability (LMI; ^3.5(1 . 
Where there is little evolution at the boundary and the ker- 
nel is chosen to be stable to transverse waves, the banding 
disappears, as was the case for our extra TSPH-HOCT4-96 
simulation. Where there is strong evolution at the bound- 
ary, as was the case for TSPH-HOCT4-442, the LMI drives 
banding irrespective of the choice of kernel. Only in our full 
scheme, OSPH-HOCT4-442, where the LMI is cured and the 
kernel is stable to transverse waves, is the banding cured. 



5.2.3 The Eo error 

Away from boundaries, the Eo error in TSPH decreases with 
the neighbour number, as expected for smooth flow (see Ap- 
pendix However, on the boundary the |Eo| error grows 
by 2-3 orders of magnitude. Increasing the neighbour num- 
ber does result in better long-term evolution, but the re- 
sults improve very slowly. This is shown in Figure [B] Notice 
that TSPH-HOCT4-442 resolves two wraps of the KH roll at 
Tkh = 2, whereas TSPH-CT-128 only manages one. How- 
ever, even in TSPH-HOCT4-442, the long-term evolution 
eventually degrades. By tkh ~ 3, the results are 'gloopy', 
rather similar to sim ulations that exp licitly model fluid sur- 
face tension (see e.g. lHerrmann|[2005^ . 

The poor Eo on the boundary is the result of a poor 
volume estimate for each particle rrij / pj (see H3.3p . However 
Eo is not solely responsible for the gloopy behaviour. There 
is a second problem - similar to a numerical surface tension 
term - that needs to be solved extra to minimising Eo . This 
is the local mixing instability (LMI) error ( i]3.5|) . 



5.2.4 The local mixing instability error 

The right panels of Figure [3] show the pressure as a function 
of y in a slice of width dx = 1 about the z-axis and width 
dx = 1 about the x-axis. In SPH and TSPH there is a clear 
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pressure discontinuity on the boundary. This is caused by 
the local mixing instability (LMI) discussed in i|3.5l 

Notice that the pressure blip is larger in TSPH than in 
SPH, yet the KHI roll progresses further in TSPH than in 
SPH. This apparent paradox is the result of the improved 
performance in TSPH. As the KHI roll progresses in TSPH, 
particles are pushed closer to the boundary making the LMI 
worse and increasing the pressure blip. In SPH, there is a 
larger gap at the boundary due to the larger surface tension 
error. This leads to less evolution and a smaller associated 
pressure blip. We will see a similar effect occurring in the 
blob test in ^ 

As discussed in H3.51 the LMI should be cured by the RT 
density estimate (equation I35|) . This is shown in Figure U 
third row which shows the results for our full OSPH scheme. 
With the RT densities, the pressure at the boundary has a 
much smaller blip, while Eq is reduced by over an order of 
magnitude. This latter effect occurs since the RT densities 
also give an improved volume estimate for each particle (see 
i|3.3l and i]3.5p . The long term evolution is now in excellent 
agreement with the RAMSES results (compare Figure [4] third 
and bottom rows). 

Although OSPH gives significantly improved results as 
compared with SPH, the scheme is numerically expensive. 
Simulations with larger density gradients require very high 
resolution. This is shown in Figure [G] second from bottom 
row. This shows the long term evolution of a KHI test with 
density ratio -Rp = 8 in OSPH. The solution should be sim- 
ilar to the Rp = 2 simulation, but it is not. The 'gloopy' 
behaviour indicative of large surface tension errors has re- 
turned. Further increasing the neighbour numbers would re- 
duce this problem, but at increased numerical cost. We will 
address this issue in future work (Hayfield & Read in prep.). 



6 THE SOD SHOCK TUBE 

Before we embark on the blob test in [JT] it is worth checking 
that our new OSPH scheme can still correctly resolve sh ocks. 
To te st this, we use a standard Sod shock tube test jSodl 

[Hzi). 

The Sod shock tube consists of a ID tube on the interval 
[—0.5, 0.5] with a discontinuous change in properties at x = 
designed to generate a shock. The left state is described by 
pi = 1.0, Pi = 1.0, vi = 0, and the right state by pr = 
0.125, Pr = 0.1, Vr = 0, where p,P and v are the density, 
pressure and velocity along the x axis. We use an adiabatic 
equation of state with 7 = 1.4. The subsequent evolution of 
the problem has a self-similar analytic solution that has a 
number distinct features which quite generally test a code's 
conservation properties, artificial viscosity, ability to handle 
nonlinear waves, and shock resolution. 

Figure [7] shows the results for the Sod shock tube test 
at time t — 0.2 in SPH (top) OSPH (bottom). Since we are 
primarily concerned with the 3D performance of the code, 
the test was performed in 3D on the union of a 24 x 24 x 300 
lattice on the left, with a 12 x 12 x 150 lattice on the right, 
giving a ID resolution of 450 points. We use 442 neighbours 
for this test in both SPH and OSPH to ensure that any 
difference is not simply due to improved kernel sampling in 
OSPH. 

For SPH, the only strong disagreement with the analytic 



solution is in the pressures that have a blip at a; = 0.2, and 
the temperatures that overshoot at 2: = 0.2. The former 
feature is due to the LMI (see 333] and i|5.2.4|) . The latter 
feature is seen in all SPH Sod shock tube tests and results 
from the well-known 'wall heating' effect (|Nohlll987l ). This 
is an error due to the artificial viscosity prescription and is 
beyond the scope of this present work. 

For OSPH, the results are even better than for SPH. 
The pressure blip is now gone, while the temperature over- 
shoot at a; = 0.2 is reduced. Only the velocities appear to be 
worse, with some remaining dispersion at x = 0.2. This owes 
to the jump in density at this point, and the associated jump 
in |Eo|. This gives a force error at the discontinuity which 
introduces some dispersion into the velocities. In SPH this 
cannot occur since the LMI causes a pressure blip at the 
boundary that prevents mixing. We will discuss this issue 
further in a forthcoming paper (Hayfield & Read in prep.). 



7 THE BLOB TEST 

The KHI test presented in ij5] is a worst-case scenario for 
OSPH, since it has a pure adiabatic sharp boundary. For 
many practical situations, boundaries will be less sharp, 
while physical entropy generation due to shocks and/or 
cooling will suppress the LMI. We give a practical exam- 
ple of this in this section using the blob test described in 
lAgertz et al.l (|2007l ') . A spherical cloud of gas of radius _Rci is 
placed in a wind tunnel with periodic boundary conditions. 
The ambient medium is ten times hotter and ten times less 
dense than the cloud so that it is in pressure equilibrium 
with the latter. We refer to the initial density contrast be- 
tween the cloud and the medium as Rp,ini- The wind velocity 
(«wind ~ CsM) has au associated Mach number M = 2.7. 
This leads to the formation of a bow shock after which the 
post-shock subsonic fiow interacts with the cloud and turns 
supersonic as it fiows past it. 

The blob test is useful for investigating how different hy- 
drodynamics codes model astrophysical processes important 
for multiphase systems, such as shocks, ram-pressure strip- 
ping and fragmentation through KH and Rayleigh- Taylor 
(RT) instabilities. As tkh < trt, an approximate timescale 
of the cloud destruction is that of the full growth of the 
largest KH mode, i.e. the wavelength of the cloud's radius. 
This can be obtained by considering the post-shock flow on 
the cloud and its time-dependence as the shock weakens and 
the cloud is accelerated. A f ull analysis of this test is pre- 
sented in lAgertz et all (|2007l ') and gives tkh ~ 1.6rcr where 
Tcr = 27?ciJ?p'ini/v is the crushing time and the velocity v 
refers to the streaming velocity in the reference frame of the 
cloud. After this time, the cloud is expected to show a more 
complicated non-linear behaviour leading to disruption. The 
original blob test was initialised in a glass-like configuration 
obtained using simulated annealing using a standard SPH 
code. Since we now use OSPH rather than SPH, we must set 
up new ICs for the blob. Our new IC set up is described in 
detail in Appendix [B] Unlike the previous blob test, where 
perturbations were seeded by random noise in the particle 
distribution, here we deliberately seed an inward growing 
mode on the front surface of the blob. This makes com- 
parison between OSPH and FLASH simpler, since then the 
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morphology of the blob is less dependent on small scale nu- 
merical noise. 

The results are presented in Figure [8l where we com- 
pare SPH, TSPH and OSPH with increasing r esolution with 
similar results from the Eulerian code FLASH (jFrvxell et alj 
|2000D . The SPH resuhs ( t op pa nels) are similar to those 
presented in lAgertz et alJ l|2007l ) . The blob is squashed by 
the shock, but does not break up. There are no visible sur- 
face KHI or Rayleigh- Taylor instabilities. TSPH (second 
row) gives significantly improved results. The central de- 
pression is now resolved and the blob is mostly destroyed 
by tkh = 3. However, the density remains clumpy as com- 
pared to the FLASH simulation (bottom row) . Our full OSPH 
scheme (third row) gives excellent agreement with the FLASH 
results. There are clear surface KHI and RT instabilities and 
the blob breaks up fully by tkh = 3. The precise details 
of the break up in FLASH and OSPH are different. How- 
ever, these differences are smaller than those we observed 
between FLASH simulations of varying resolution. They are 
caused by the non-linear break up of the blob that is affected 
by resolution-scale perturbations. 

Figure |5] shows Eo and the pressure blips for the blob 
test in SPH (left), TSPH (middle) and OSPH (right) at 
Tkh = 1. In SPH and TSPH, the two fluid phases (marked 
by the black and grey solid circles) remain well separated 
at all times. In both cases the pressure distribution shows 
discontinuities. By contrast, in OSPH the fluids are already 
mixed at tkh = 1, while the pressures are smooth and single- 
valued throughout the flow. 

There is a more modest improvement in Eq between 
SPH, TSPH and OSPH than that seen in the KHI tests 
presented in ^ OSPH gives a Eo smaller by a factor ~ 5 
compared to the SPH and TSPH simulations, whereas in 
the KHI tests, there was an improvement of over an order of 
magnitude. There are two reasons for this. Firstly, since the 
SPH simulation shows little evolution at the boundary, the 
initial Eq is relatively well conserved. By contrast, TSPH 
shows significant boundary evolution due to improved mix- 
ing. This can actually worsen Eq at the boundary since the 
particles are still unable to properly interpenetrate as a re- 
sult of the LMI. Figure [8] clearly shows, however, that TSPH 
gives improved mixing. Secondly, entropy generation at the 
shock softens the density step around the blob, leading to 
an improved Eq even for the SPH simulation. 



8 CONCLUSIONS 

Standard formulations of smoothed particle hydrodynamics 
(SPH) cannot resolve fluid mixing and instabilities at flow 
boundaries. We have used an error and stability analysis of 
the generalised SPH equations of motion to show that mix- 
ing fails for two distinct reasons. The first is a leading order 
error in the momentum equation. This should decrease with 
increasing neighbour number, but does not because numer- 
ical instabilities cause the kernel to be irregularly sampled. 
We identified two important instabilities: the clumping in- 
stability and the banding instability, and we showed that 
both are cured by a suitable choice of kernel. The second 
problem is the local mixing instability (LMI). This occurs 
as particles attempt to mix on the kernel scale, but are un- 
able to due to entropy conservation. The result is a pressure 



discontinuity at boundaries that pushes fluids of different en- 
tropy apart. We cure d the LMI by using a weig hted density 
estimate proposed bv lRitchie fc Thomaj l|2001 ). We showed 
that this both reduces errors in the continuity equation and 
allows individual particles to mix at constant pressure. 

We demonstrated mixing in our new Optimised 
Smoothed Particle Hydrodynamics (OSPH) scheme using a 
Kelvin Helmholtz instability (KHI) test with density con- 
trast 1:2, and the 'blob test' - a 1:10 density ratio gas 
sphere in a wind tunnel - finding excellent agreement be- 
tween OSPH and Eulerian codes. 

OSPH is a multiphase Lagrangian method that con- 
serves momentum, mass and entropy, and demonstrates that 
it is possible to model multiphase fluid flow using SPH. How- 
ever, OSPH remains a low-order method, requiring large 
neighbour number to keep the Eo error small. We will ad- 
dress this problem in a forthcoming paper, where we use the 
lessons learnt in this present work to move to higher order 
particle methods (Hayfield & Read in prep.). 
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APPENDIX A: THE EFFECT OF INCREASING 
NEIGHBOUR NUMBER IN TSPH 

In this appendix we show the effect of increasing neighbour 
number for TSPH-CT (i.e. without the clumping instabil- 
ity). The results for the same KH instability test shown in 
Figure |4] are shown in Figure [KT\ for 32 and 64 neighbours. 
Notice that with increasing particle number, the error vector 
Eo is reduced, the pressure blip at the boundary is reduced, 
and the results improve. 



APPENDIX B: THE BLOB TEST SETUP 

The hydrodynamical properties of the blob test are de- 
scribed in SJT] We use a periodic simulation box of size, in 
units of the cloud radius _Rci, {L^, Ly, L^} = {10,10,30} 
and we centre the cloud at {x,y,z} = {5,5,5}. The equal 
mass SPH particles constituting the ambient medium and 
the cloud are arranged in lattice configurations to achieve 
the relevant density contrast i?p, ini ~ 10. The particle tem- 
peratures (T ~ P/p) are then assigned to achieve pressure 
equilibrium where the local density measurement of equation 
[T]is used for consistency. The wind velocity (wwind = CsM) 
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TSPH-CT-32 




Figure Al. As FigureU but showing the effect of increasing particle number for TSPH-CT. 



has an associat ed Ma ch number M — 2.7, where the sound 
speed is Cs = \J^PI p using an adiabatic index 7 = 5/3. 

In the original blob test described in lAgertz et al.l 
SPH particle noise was used to trigger instabili- 
ties. This procedure is not applicable when using a noise 
free lattice configuration. Hence, we use spherical harmon- 
ics to apply large scale perturbations to the surface layer 
of the cloud. The full perturbation, in spherical coordinates 
centred on the cloud, can be expressed as Upcrt(r, 6, 0) = 
5«i?(r)Re[y(6l,<^)r)]/C, where the radial part, R{r) = 
exp (2(r — rci)/rci), is defined for r ^ rd and the spherical 
harmonic is, adopting Z = 5, m = 3: 

32 y TT 

The constant C simply normalises the real part of the har- 
monic to reach a maximum value of 1. We chose a subsonic 
perturbation &v — — 0.06uwind. 
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